In a slight detour in this project, this section concerns fitting time series data with statistical models that concern the variance in the data. These models are called AutoRegressive Conditional Heteroskedasticity (ARCH) and Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) models and are used when there is an autoregressive term in the error variance. They are typically used when modeling time-series data that exhibit volatility, which happens most often in financial time series. Therefore, for this section, a set of stock data is retrieved and analyzed using a combination of ARCH/GARCH and ARIMA.
NextEra Energy Inc
For this section, the stock of choice is NextEra Energy Inc, which is quite volatile and connected to the renewable energy as it has the largest market cap or net worth of 159.73 billion U.S. dollars among other renewable energy companies. After gathering the adjusted closing prices from Yahoo Finance, the following visualization below shows the price trend over 20+ years.
Stock Prices
Gather NextEra Energy stocks from Yahoo Fiance
NEE is the ticker symbol for NextEra Energy’s stock.
Code
NEE_ALL <-getSymbols("NEE",auto.assign =FALSE, from ="2019-01-01",to ="2023-04-13",src="yahoo")NEE_ALL=data.frame(NEE_ALL)NEE_ALL <-data.frame(NEE_ALL, rownames(NEE_ALL))colnames(NEE_ALL)[7] ="Year"NEE_ALL$Year <-as.Date(NEE_ALL$Year,"%Y-%m-%d")knitr::kable(head(NEE_ALL,3))
# plotly# change font family and sizet <-list(family ="Palatino",size =14)theme_set(theme_gray(base_size=12,base_family="Palatino"))fig <- NEE_ALL |>plot_ly(x =~Year, type="candlestick",open =~NEE.Open, close =~NEE.Close,high =~NEE.High, low =~NEE.Low, name="NEE") fig <- fig |>add_lines(x =~Year, y =ma(NEE_ALL$NEE.Adjusted, 20), name="20-MA", line =list(color ='deeppink', width =1), inherit = F)fig <- fig |>add_lines(x =~Year, y =ma(NEE_ALL$NEE.Adjusted, 80), name="80-MA", line =list(color ='cyan3', width =1), inherit = F)fig <- fig |>layout(title ="Candlestick Chart of NextEra Energy Stock Prices",font=t, legend =list(orientation ='h', x =0.5, y =1,xanchor ='center', yref ='paper', font =list(size =10), bgcolor ='transparent'),xaxis =list(type ='date',tickformat ="%d %B (%a)<br>%Y" ),plot_bgcolor ="#eeeeee" )fig
Code
p <- NEE_ALL |>ggplot() +geom_line(aes(y=NEE.Adjusted, x=Year),color="aquamarine3") +ggtitle("NextEra Energy Stock Pirces") +ylab("USD") +xlab("Month") +scale_x_date(date_breaks ="1 month",labels = scales::date_format(format ="%Y-%m")) +# rotate labels for visibilitytheme(axis.text.x =element_text(angle =270, vjust =0.5, hjust=1))ggplotly(p)
From the above figures, we can see the spikes of volatility in the price, especially around certain points where the price dramatically drops and rises (see around mid 2020). This is a good first indicator that there will be an ARCH/GARCH component when fitting the data. Another way to look at this data is by taking a look at price versus volume. This just gives us a better sense of the price and a different perspective, but it doesn’t necessarily provide us with more insight.
Price Volumn Chart
Code
NEE_ALL <-getSymbols("NEE",auto.assign =FALSE, from ="2019-01-01",to ="2023-04-13",src="yahoo")chartSeries(NEE_ALL, theme =chartTheme("white"), # Themebar.type ="hlc", # High low close up.col ="green", # Up candle colordn.col ="red",name ="NEE Stock") # Down candle color)
Again, the above figure shows us more information on the price trends versus the volume, but the overall conclusion is the same. However, the next figure shows the volatility better as the log of the returns is taken and plotted.
Code
returns <-diff(log(NEE_ALL$`NEE.Adjusted`), lag =1)returns %>%chartSeries(theme =chartTheme("white"), name ="NEE Stock Daily Log Returns")
The above figure shows the daily log of the returns, which gives a much better sense of the volatility of the returns for NextEra Energy. Here, since there is clear volatility, an ARCH/GARCH model will be needed to fit the data. However, before fitting an ARCH model immediately, it is important to check the ACF/PACF plots to check if an ARIMA model should be fit first before fitting an ARCH model on the residuals. First, since we took the log of the returns, we should plot the ACF to check if we should difference the data.
ACF and PACF Plots
Code
log(NEE_ALL$'NEE.Adjusted') |>ggAcf() +ggtitle("ACF of Log Transformed NextEra Energy Prices")
From the above figure, it is clear that a difference must be done to account for the correlation in the data. After differencing the data once, the ACF and PACF plots are created and shown below.
acf <-ggAcf(returns, lag.max =100) +ggtitle("ACF of Log Transformed and Differenced NextEra Energy Prices")ggplotly(acf)
Code
pacf <-ggPacf(returns, lag.max =100) +ggtitle("PACF of Log Transformed and Differenced NextEra Energy Prices")ggplotly(pacf)
Seasonality Test
Code
isSeasonal(returns, test ="combined", freq =NA)
[1] FALSE
According to the seasonality test, NextEra stock price is not seasonal time series data.
Model Fitting - ARIMA Model
Now that the data is differenced, the ACF and PACF plots clearly show a need for an ARIMA model to fit the data. Once an ARIMA model is fit, there will probably be a need to fit an ARCH model on the residuals, but that will depend on the residuals. Let’s fit an ARIMA model first. From the above ACF and PACF plots, it would be appropriate to check all Arima models with p values ranging from 1-7 and q values ranging from 1-7. For the ARIMA model, d is set to 1 as we differenced the data once. From that, the following results were achieved.
Code
d =1i =1temp =data.frame()ls =matrix(rep(NA,6*43),nrow=43)for (p in1:7){for(q in1:7) {if(p + d + q <=12) { model <-Arima(log(NEE_ALL$'NEE.Adjusted'),order=c(p,d,q), include.drift =TRUE) ls[i,] =c(p,d,q,model$aic,model$bic,model$aicc) i = i +1 } }}temp =as.data.frame(ls)names(temp) =c("p","d","q","AIC","BIC","AICc")
By sorting the table, the smallest AIC has the parameters of 4,1,5, while by sorting using the BIC criterion, the model has parameters of 2,1,2. Next, a function allows for an auto-selection of an ARIMA model. According to the auto.arima() method in R, the parameters for the model should be 2,1,2. Therefore, the 2,1,2 model makes the most sense.
Next, the model should undergo some diagnostics to make sure the fit looks good. Below are some basic diagnostics run using the 2,1,2 model.
Given the figure above, the following observations can be made. First, when looking at the plot of standardized residuals, the mean should be around 0, with the variance at about 1. As seen above, this is generally true, with the mean close to 0. However, the variance is probably slightly more than 1 and there is clusters of higher variance, which indicates a second model should be fit for the errors. Generally, if the mean or variance is significantly different than expected, this would be a sign of a poor model fit, but we don’t really see that here. Next, when inspecting the ACF of residuals, the plot shows no significant lags, which is a very encouraging sign considering the logic above. Next, when looking at the qq-plot, the hope is to see some signs of normality. Here, we see some signs of normality in the data, but with a bit of skew. Finally, when looking at the p-values for the Ljung-Box test, many of the p-values are above 0.05, which is a decent sign of a good model.
Next, we can do an initial check on the absolute returns and squared returns to check if an ARCH model is needed. Below is a set of visualizations for the ACF and PACF plots for the absolute and squared returns.
fit <-arima(log(NEE_ALL$'NEE.Adjusted'), order =c(2, 1, 2))res.arima <- fit$ressquared.res.arima <- res.arima^2acf_residuals <-ggAcf(squared.res.arima, na.action = na.pass) +ggtitle("ACF Squared Residuals")ggplotly(acf_residuals)
Code
fit <-arima(log(NEE_ALL$'NEE.Adjusted'), order =c(2, 1, 2))res.arima <- fit$ressquared.res.arima <- res.arima^2pacf_residuals <-ggPacf(squared.res.arima, na.action = na.pass) +ggtitle("ACF Squared Residuals")ggplotly(pacf_residuals)
From the above plots and the previous diagnostics, it is pretty clear an ARCH model is needed to fit the residuals of the ARIMA model. Using the above plots, the p-value should range between 1-9. Then, using a similar technique to fitting the ARIMA model, a set of ARCH models can be created and evaluated to find the best fit for the residuals. Below is another table that shows the fitted ARCH models with their associated AICs.
Model Fitting - ARIMA(2,1,2) + ARCH Model
Code
ARCH <-list() ## set countercc <-1for (p in1:9) { ARCH[[cc]] <-garch(res.arima, order =c(0, p), trace = F) cc <- cc +1}## get AIC values for model evaluationARCH_AIC <-sapply(ARCH, AIC) ## model with lowest AIC is the bestmin(ARCH_AIC)
By sorting the table, the smallest AIC has the parameters of 6,0. Therefore, we should fit the residuals with a 6,0 ARCH model. Now that we have fit the data with an ARCH/ARIMA model, the next step is to check the fit. One option is to check the residuals. First, we can plot the residuals and evaluate the results. Additionally, we can check whether the residuals are normal.
From above, we can see a couple of important diagnostic plots. The first set of plots show a set of 3 plots of the residuals, the ACF of the residuals, and a distribution of the residuals. The residuals seem to now have a constant mean and variance according to the first plot. Next, the ACF shows strong signs of stationarity in the residuals. Finally, the distribution of the residuals is relatively normal. This can be confirmed using the second tab, which is a qqnorm plot of the residuals. Using this plot, the normality of the residuals can be confirmed. These plots all point towards strong signs the model is a good fit for the data. Finally, a Ljung-Box test outputs a p-value greater than 0.05, indicating that the residuals are independent, encouraging the suitable model. Here, after doing the model diagnostics, the model above can be used for forecasting, and the model is ARCH(6,0) + ARIMA(2,1,2).
Using this model, the data can be effectively modeled and perhaps future values can be forecasted. However, stock market prices are generally impossible to forecast. We can say this model does a decent job at modeling the data and its variance.
Source Code
---title: "Financial Time Series Models (ARCH/GARCH)"format: html: embed-resources: true self-contained: true code-fold: true code-copy: true code-tools: true code-overflow: wrap---```{r,echo=FALSE,message=FALSE,warning=FALSE}library(tidyverse)library(ggplot2)library(forecast)library(astsa) library(xts)library(tseries)library(fpp2)library(fma)library(lubridate)library(tidyverse)library(TSstudio)library(quantmod)library(tidyquant)library(plotly)library(ggplot2)library(ggfortify)library(autoplotly)library(ggpubr)library(seastests)library(fGarch) library(dynlm)theme_set(theme_gray(base_size=12,base_family="Palatino"))```In a slight detour in this project, this section concerns fitting time series data with statistical models that concern the variance in the data. These models are called AutoRegressive Conditional Heteroskedasticity (ARCH) and Generalized AutoRegressive Conditional Heteroskedasticity (GARCH) models and are used when there is an autoregressive term in the error variance. They are typically used when modeling time-series data that exhibit **volatility**, which happens most often in financial time series. Therefore, for this section, a set of stock data is retrieved and analyzed using a combination of ARCH/GARCH and ARIMA.# NextEra Energy Inc{fig-align="center" width="525"}For this section, the stock of choice is [**NextEra Energy Inc**](https://www.google.com/finance/quote/NEE:NYSE?sa=X&ved=2ahUKEwiGuv6Whbj-AhVyk2oFHUk-ATQQ3ecFegQILBAg), which is quite volatile and connected to the renewable energy as it has the largest market cap or net worth of 159.73 billion U.S. dollars among other renewable energy companies. After gathering the adjusted closing prices from Yahoo Finance, the following visualization below shows the price trend over 20+ years.## Stock Prices## Gather NextEra Energy stocks from Yahoo Fiance`NEE` is the ticker symbol for NextEra Energy's stock.```{r}NEE_ALL <-getSymbols("NEE",auto.assign =FALSE, from ="2019-01-01",to ="2023-04-13",src="yahoo")NEE_ALL=data.frame(NEE_ALL)NEE_ALL <-data.frame(NEE_ALL, rownames(NEE_ALL))colnames(NEE_ALL)[7] ="Year"NEE_ALL$Year <-as.Date(NEE_ALL$Year,"%Y-%m-%d")knitr::kable(head(NEE_ALL,3))```## Visualize NEE Stock::: panel-tabset## Candlestick Plot```{r}# plotly# change font family and sizet <-list(family ="Palatino",size =14)theme_set(theme_gray(base_size=12,base_family="Palatino"))fig <- NEE_ALL |>plot_ly(x =~Year, type="candlestick",open =~NEE.Open, close =~NEE.Close,high =~NEE.High, low =~NEE.Low, name="NEE") fig <- fig |>add_lines(x =~Year, y =ma(NEE_ALL$NEE.Adjusted, 20), name="20-MA", line =list(color ='deeppink', width =1), inherit = F)fig <- fig |>add_lines(x =~Year, y =ma(NEE_ALL$NEE.Adjusted, 80), name="80-MA", line =list(color ='cyan3', width =1), inherit = F)fig <- fig |>layout(title ="Candlestick Chart of NextEra Energy Stock Prices",font=t, legend =list(orientation ='h', x =0.5, y =1,xanchor ='center', yref ='paper', font =list(size =10), bgcolor ='transparent'),xaxis =list(type ='date',tickformat ="%d %B (%a)<br>%Y" ),plot_bgcolor ="#eeeeee" )fig```## Time Series Plot```{r}p <- NEE_ALL |>ggplot() +geom_line(aes(y=NEE.Adjusted, x=Year),color="aquamarine3") +ggtitle("NextEra Energy Stock Pirces") +ylab("USD") +xlab("Month") +scale_x_date(date_breaks ="1 month",labels = scales::date_format(format ="%Y-%m")) +# rotate labels for visibilitytheme(axis.text.x =element_text(angle =270, vjust =0.5, hjust=1))ggplotly(p)```:::From the above figures, we can see the spikes of volatility in the price, especially around certain points where the price dramatically drops and rises (see around mid 2020). This is a good first indicator that there will be an ARCH/GARCH component when fitting the data. Another way to look at this data is by taking a look at price versus volume. This just gives us a better sense of the price and a different perspective, but it doesn't necessarily provide us with more insight.### Price Volumn Chart```{r}NEE_ALL <-getSymbols("NEE",auto.assign =FALSE, from ="2019-01-01",to ="2023-04-13",src="yahoo")chartSeries(NEE_ALL, theme =chartTheme("white"), # Themebar.type ="hlc", # High low close up.col ="green", # Up candle colordn.col ="red",name ="NEE Stock") # Down candle color)```Again, the above figure shows us more information on the price trends versus the volume, but the overall conclusion is the same. However, the next figure shows the volatility better as the log of the returns is taken and plotted.```{r}returns <-diff(log(NEE_ALL$`NEE.Adjusted`), lag =1)returns %>%chartSeries(theme =chartTheme("white"), name ="NEE Stock Daily Log Returns")```The above figure shows the daily log of the returns, which gives a much better sense of the volatility of the returns for NextEra Energy. Here, since there is clear volatility, an ARCH/GARCH model will be needed to fit the data. However, before fitting an ARCH model immediately, it is important to check the ACF/PACF plots to check if an ARIMA model should be fit first before fitting an ARCH model on the residuals. First, since we took the log of the returns, we should plot the ACF to check if we should difference the data.## ACF and PACF Plots```{r}log(NEE_ALL$'NEE.Adjusted') |>ggAcf() +ggtitle("ACF of Log Transformed NextEra Energy Prices")```From the above figure, it is clear that a difference must be done to account for the correlation in the data. After differencing the data once, the ACF and PACF plots are created and shown below.::: panel-tabset## ACF```{r}acf <-ggAcf(returns, lag.max =100) +ggtitle("ACF of Log Transformed and Differenced NextEra Energy Prices")ggplotly(acf)```## PACF```{r}pacf <-ggPacf(returns, lag.max =100) +ggtitle("PACF of Log Transformed and Differenced NextEra Energy Prices")ggplotly(pacf)```:::### Seasonality Test```{r}isSeasonal(returns, test ="combined", freq =NA)```According to the seasonality test, NextEra stock price is not seasonal time series data.## Model Fitting - ARIMA ModelNow that the data is differenced, the ACF and PACF plots clearly show a need for an ARIMA model to fit the data. Once an ARIMA model is fit, there will probably be a need to fit an ARCH model on the residuals, but that will depend on the residuals. Let's fit an ARIMA model first. From the above ACF and PACF plots, it would be appropriate to check all Arima models with p values ranging from 1-7 and q values ranging from 1-7. For the ARIMA model, d is set to 1 as we differenced the data once. From that, the following results were achieved.```{r}d =1i =1temp =data.frame()ls =matrix(rep(NA,6*43),nrow=43)for (p in1:7){for(q in1:7) {if(p + d + q <=12) { model <-Arima(log(NEE_ALL$'NEE.Adjusted'),order=c(p,d,q), include.drift =TRUE) ls[i,] =c(p,d,q,model$aic,model$bic,model$aicc) i = i +1 } }}temp =as.data.frame(ls)names(temp) =c("p","d","q","AIC","BIC","AICc")``````{r}#knitr::kable(temp)knitr::kable(temp[which.min(temp$AIC),])knitr::kable(temp[which.min(temp$BIC),])# auto.arima()auto.arima(log(NEE_ALL$'NEE.Adjusted'))```By sorting the table, the smallest AIC has the parameters of 4,1,5, while by sorting using the BIC criterion, the model has parameters of 2,1,2. Next, a function allows for an auto-selection of an ARIMA model. According to the `auto.arima()` method in R, the parameters for the model should be 2,1,2. Therefore, the 2,1,2 model makes the most sense.Next, the model should undergo some diagnostics to make sure the fit looks good. Below are some basic diagnostics run using the 2,1,2 model.```{r, warning=FALSE}sarimafit <-sarima(log(NEE_ALL$'NEE.Adjusted'), 2, 1, 2, details =FALSE)```{fig-align="center"}Given the figure above, the following observations can be made. First, when looking at the plot of standardized residuals, the mean should be around 0, with the variance at about 1. As seen above, this is generally true, with the mean close to 0. However, the variance is probably slightly more than 1 and there is clusters of higher variance, which indicates a second model should be fit for the errors. Generally, if the mean or variance is significantly different than expected, this would be a sign of a poor model fit, but we don't really see that here. Next, when inspecting the ACF of residuals, the plot shows no significant lags, which is a very encouraging sign considering the logic above. Next, when looking at the qq-plot, the hope is to see some signs of normality. Here, we see some signs of normality in the data, but with a bit of skew. Finally, when looking at the p-values for the Ljung-Box test, many of the p-values are above 0.05, which is a decent sign of a good model.Next, we can do an initial check on the absolute returns and squared returns to check if an ARCH model is needed. Below is a set of visualizations for the ACF and PACF plots for the absolute and squared returns.::: panel-tabset## Squared Residuals ACF```{r}fit <-arima(log(NEE_ALL$'NEE.Adjusted'), order =c(2, 1, 2))res.arima <- fit$ressquared.res.arima <- res.arima^2acf_residuals <-ggAcf(squared.res.arima, na.action = na.pass) +ggtitle("ACF Squared Residuals")ggplotly(acf_residuals)```## Squared Residuals PACF```{r}fit <-arima(log(NEE_ALL$'NEE.Adjusted'), order =c(2, 1, 2))res.arima <- fit$ressquared.res.arima <- res.arima^2pacf_residuals <-ggPacf(squared.res.arima, na.action = na.pass) +ggtitle("ACF Squared Residuals")ggplotly(pacf_residuals)```:::From the above plots and the previous diagnostics, it is pretty clear an ARCH model is needed to fit the residuals of the ARIMA model. Using the above plots, the p-value should range between 1-9. Then, using a similar technique to fitting the ARIMA model, a set of ARCH models can be created and evaluated to find the best fit for the residuals. Below is another table that shows the fitted ARCH models with their associated AICs.## Model Fitting - ARIMA(2,1,2) + ARCH Model```{r}ARCH <-list() ## set countercc <-1for (p in1:9) { ARCH[[cc]] <-garch(res.arima, order =c(0, p), trace = F) cc <- cc +1}## get AIC values for model evaluationARCH_AIC <-sapply(ARCH, AIC) ## model with lowest AIC is the bestmin(ARCH_AIC)which(ARCH_AIC ==min(ARCH_AIC))ARCH[[which(ARCH_AIC ==min(ARCH_AIC))]]arch.df <-data.frame(p =1:9, q =0, AIC = ARCH_AIC)knitr::kable(arch.df)```By sorting the table, the smallest AIC has the parameters of 6,0. Therefore, we should fit the residuals with a 6,0 ARCH model. Now that we have fit the data with an ARCH/ARIMA model, the next step is to check the fit. One option is to check the residuals. First, we can plot the residuals and evaluate the results. Additionally, we can check whether the residuals are normal.::: panel-tabset## Diagnostics```{r, warning=FALSE}fit2 <-garch(res.arima, order =c(0,6), trace =FALSE)checkresiduals(fit2)```## Normality Tests```{r}qqnorm(fit2$residuals, pch =1)qqline(fit2$residuals, col ="blue", lwd =2)```## Ljung-Box Test```{r}Box.test(fit2$residuals, type ="Ljung")```:::From above, we can see a couple of important diagnostic plots. The first set of plots show a set of 3 plots of the residuals, the ACF of the residuals, and a distribution of the residuals. The residuals seem to now have a constant mean and variance according to the first plot. Next, the ACF shows strong signs of stationarity in the residuals. Finally, the distribution of the residuals is relatively normal. This can be confirmed using the second tab, which is a qqnorm plot of the residuals. Using this plot, the normality of the residuals can be confirmed. These plots all point towards strong signs the model is a good fit for the data. Finally, a Ljung-Box test outputs a p-value greater than 0.05, indicating that the residuals are independent, encouraging the suitable model. Here, after doing the model diagnostics, the model above can be used for forecasting, and the model is **ARCH(6,0) + ARIMA(2,1,2)**.## Forecast```{r, warning=FALSE}fit3 =garchFit(~arma(2,2) +garch(6, 0), data = res.arima, trace =FALSE)predict(fit3, n.ahead =5, plot=TRUE)```Using this model, the data can be effectively modeled and perhaps future values can be forecasted. However, stock market prices are generally impossible to forecast. We can say this model does a decent job at modeling the data and its variance.
